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Meron-Cluster Solution of Fermion and Other Sign Problems 
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Numerical simulations of numerous quantum systems suffer from the notorious sign problem. Important ex- 
amples include QCD and other field theories at non-zero chemical potential, at non-zero vacuum angle, or with 
an odd number of flavors, as well as the Hubbard model for high-temperature superconductivity and quantum 
antiferromagnets in an external magnetic field. In all these cases standard simulation algorithms require an ex- 
ponentially large statistics in large space-time volumes and are thus impossible to use in practice. Meron-cluster 
algorithms realize a general strategy to solve severe sign problems but must be constructed for each individual 
case. They lead to a complete solution of the sign problem in several of the above cases. 
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1. Introduction 

Sign problems prevent the numerical solution 
of several important problems in physics. For ex- 
ample, QCD and other field theories at non-zero 
chemical potential or at non-zero vacuum angle 
have a complex action and hence complex Boltz- 
mann weights which cannot be interpreted as 
probabilities in a Monte Carlo calculation. When 
the complex phase of the Boltzmann factor is in- 
cluded in measured observables, the fluctuations 
of the phase give rise to dramatic cancellations. 
Especially for large systems at low temperatures 
this leads to relative statistical errors that are 
exponentially large in both the volume and the 
inverse temperature. Similarly, the minus-signs 
that arise as a consequence of Fermi statistics 
cause severe sign problems in numerical simula- 
tions of strongly correlated electron systems such 
as the Hubbard model for high-temperature su- 
perconductors. Another class of systems that suf- 
fer from a sign problem are frustrated quantum 
spin systems — for example, quantum antiferro- 
magnets in an external magnetic field. In all these 
cases it is impossible in practice to study these 
systems with standard numerical methods. 

Recently, it has been shown to be possible to 
completely solve some severe sign problems us- 
ing meron-cluster algorithms [0J|]. The solution 
of the problem proceeds in two steps. The idea 
of the first step is to decompose a field config- 



uration into independent clusters whose fiip af- 
fects the sign in a definite way. Clusters whose 
flip changes the sign are referred to as merons. 
The flip of a meron-cluster leads to an exact can- 
cellation between two contributions ±1. Using 
an improved estimator, this reduces the problem 
of canceling many contributions ±1 to the prob- 
lem of averaging over non-negative contributions 
and 1. An algorithm based on this first step 
has been successfully applied to the simulation of 
a bosonic model with a complex action — the 2-d 
0(3) model at non- vanishing 0- vacuum angle d). 
In this model the sign-changing clusters are half- 
instantons — thus the name moron. However, the 
improved estimator alone solves only one half of 
the sign problem because most of the time one 
generates contributions of to the average sign 
and very rarely one encounters a contribution of 
1. In order to solve the other half of the problem 
a second step is necessary which guarantees that 
contributions of and 1 are generated with simi- 
lar probabilities. The idea behind the second step 
is to include a Metropolis decision in the process 
of cluster decomposition. The two basic ideas be- 
hind the algorithm are general and can be applied 
to a variety of systems. However, even though the 
meron concept is universal, the meron-cluster al- 
gorithm must be constructed for each individual 
case. Here we present several cases, from both 
particle and condensed matter physics, in which 
a meron-cluster algorithm has led to a complete 
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solution of a severe sign problem. 

2. General Nature of the Sign Problem 

First of all. it should be stressed that the sign 
problem is strongly influenced by the choice of ba- 
sis for the physical Hilbert space. When one con- 
structs a path integral representation of a quan- 
tum statistical partition function 



Z = Trexp(-/3if), 



(1) 



one divides the Euclidean time interval of ex- 
tent (3 into M small time-steps of size e (with 
(3 = Me) and one inserts complete sets of ba- 
sis states |n) between the operators exp(— eTJ). 
The product of the resulting transfer matrix el- 
ements (n| exp(— eiJ)|n') defines the Boltzmann 
factor. Depending on the choice of basis the 
Boltzmann weight may be positive, negative or 
even complex. When the Boltzmann weight is 
complex, one can always restrict oneself to the 
real part, because the total partition function is 
always real (as long as H is Hermitean). Still, 
the Boltzmann weight may be negative and in 
general it takes the form Sign[n] exp(— 5'[n]), with 
Sign[n] = ±1 and exp(— ^[n]) > 0. In principle, 
it is always possible to avoid the sign problem 
(i.e. one can always ensure that Sign[n] = 1) by 
a clever choice of basis |n) of the physical Hilbert 
space. For example, when one chooses the ba- 
sis of Hamiltonian eigenstates, H\n) = En\n), 
all transfer matrix elements (n| cxp(— eiJ)|n') = 
exp(— eE'rOf^run' and hence the Boltzmann weights 
are non-negative. Of course, this is a rather aca- 
demic solution of the sign problem, because the 
cases we want to simulate are the ones for which 
we don't know how to diagonalize the Hamilto- 
nian. Still, the; argument shows that a change of 
basis in the path integral may have a large impact 
on the sign problem. 

Let us now consider a general path integral 



Z = ^Sign[n]exp(-5[n]) 



(2) 



over configurations n with a Boltzmann weight of 
Sign[n] = ±1 and magnitude exp(— S'[n]). Here 
S[n] is the action of a modified model, with par- 
tition function Z' = X]„ exp(— S'[n]), which does 



not suffer from the sign problem and which can 
thus be simulated with standard Monte Carlo 
methods. An observable 0\n] of the original 
model is obtained in a simulation of the modi- 
fied ensemble as 

(O) = i^O[n]Sign[n]exp(-5[n]) = 

(3) 

The average sign in the modified ensemble is given 

by 

(Sign) = i;^Sign[n]exp(-5[n]) 



- = exp(-/3yA/). 



(4) 



The last equality points to the heart of the sign 
problem. The expectation value of the sign is ex- 
ponentially small in both the volume V and the 
inverse temperature /3 because the difference be- 
tween the free energy densities A/ = / — /' of 
the original and the modified systems is necessar- 
ily positive. 

Even in an ideal simulation of the modified 

ensemble which generates N completely uncorre- 
lated configurations, the relative statistical error 
of the sign is 



ASign _ V (Sign ) - (Sign)2 ^ exp(/3yA/) 



7V(Sign) 



N 



(5) 



Here we have used that Sign =1. To de- 
termine the average sign with sufficient accu- 
racy one needs to generate on the order of A'' = 
exp(2/3VA/) configurations. For large volumes 
and small temperatures this is impossible in prac- 
tice. It is possible to solve one half of the prob- 
lem if one can match any contribution —1 with 
another contribution 1 to give 0, such that only 
a few unmatched contributions 1 remain. Then 
effectively Sign = 0, 1 and hence Sign^ = Sign. 
This reduces the relative error to 



ASign _ y^(Sign) - (Sign)^ _ exp( 



2 



(Sign) 



/]V'(Sign) 



IN' 



(6) 



One gains an exponential factor in statistics, 
but one still needs to generate N' = VN = 



3 



exp(/3yA/) independent configurations in order 
to determine the average sign accurately. This 
difficulty still arises because one generates expo- 
nentially many vanishing contributions before one 
encounters a contribution 1. As explained below, 
in the meron-clustcr algorithm an explicit match- 
ing of contributions —1 and 1 is achieved using 
an improved estimator. This step solves one half 
of the sign problem. In a second step involving 
a Metropolis decision, the algorithm ensures that 
contributions and 1 occur with similar probabil- 
ities. This step saves another exponential factor 
in statistics and solves the other half of the sign 
problem. 

3. Chiral Phase Transition with Staggered 
Fermions ^ 

To illustrate the power of the meron-cluster 
algorithm, we apply it to a model of (3 + l)-d 
staggered fermions in the Hamiltonian formula- 
tion. The model has N = 2 flavors and a Z(2) 
chiral symmetry that is spontaneously broken at 
low temperatures. The fermion determinant can 
be negative in this model. Hence, due to the 
sign problem, standard fermion algorithms fail 
in this case. The meron-cluster algorithm is the 
only numerical method available to simulate this 
model. In our method we do not integrate out 
the fermions but describe them in a Fock state 
basis. The resulting bosonic model of fermion 
occupation numbers interacts locally, but has a 
non-local fermion permutation sign resulting from 
the Pauli exclusion principle. Standard numerical 
methods would suffer from severe cancellations of 
positive and negative contributions to the parti- 
tion function. Like other cluster methods, the 
meron-clustcr algorithm decomposes a configura- 
tion of fermion occupation numbers into clusters 
which can be flipped independently. Under a clus- 
ter flip an occupied site becomes empty and vice 
versa. The main idea of the meron-cluster algo- 
rithm is to construct the clusters such that they 
affect the fermion sign independently when they 
are flipped. In addition, it must always be pos- 
sible to flip the clusters into a reference config- 
uration with a positive sign. Like other cluster 

^Based on a talk presented by K. Holland 



algorithms, the meron algorithm substantially re- 
duces critical slowing down. This advantage al- 
lows us to work directly in the chiral limit. 

We consider staggered fermions hopping on a 
3-d cubic spatial lattice with V = sites x [L 
even) and with periodic or antiperiodic spatial 
boundary conditions. We start in the Hamilto- 
nian formulation and then derive a path integral 
on a (3 + l)-d Euclidean space-time lattice. The 
fermions are described by creation and annihila- 
tion operators and with standard anti- 
commutation relations 

{*J , *+} = * J = 0, * J = 5xy. (7) 

The staggered fermion Hamilton operator takes 

the form 

x,i X 

that is a sum of nearest-neighbor couplings hx,i 
and a mass term to\E'\1/. From here on, we work 
directly in the chiral limit m = 0, and only use 
as an observable. The term hx,i couples the 
fermion operators at the lattice sites x and x + i, 
where « is a unit- vector in the i-direction, and 

+ (9) 

Here ry^^.i = 1, ??a;,2 = (-1)''' and r/a;,3 = 
(—1)^1+^^ are the standard staggered fermion 
sign factors and G is a four-fermion coupling con- 
stant. 

To construct a path integral for the partition 
function, we decompose the Hamilton operator 
into six terms H = Hi + H2 + ... + Hq with 

Hi = ^ hx,i, = ^ hx,i- (10) 

Xieven Xiodd 

The individual contributions to a given iJj com- 
mute with each other, but two different Hi do not 
commute. Using the Suzuki- Trotter formula, we 
express the fermionic partition function at inverse 

temperature /? as 

Z = Trexp(-/3iJ) 

= lim Tr[exp(-effi)...exp(-eff6)]^-^- (H) 

M— >oo 
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We have introduced 6M Euclidean time slices 
with e — (3/M being the lattice spacing in the 
Euclidean time direction. Following Jordan and 
Wigner Q we represent the fermion operators by 
Pauli matrices 

^+ ^afc7l..c7f_,a+, C7fal..af_ia-, 
n, = vl/+vl,^ = i(af + 1), (12) 
with 

'^^ = li'^l ± ^'^f )> i'^K] = 2iSi„,e,jk<yi- (13) 

Here / labels the lattice point x. 

We now insert complete sets of fermion Fock 
states between the factors exp(— eiJj). Each site 
is either empty or occupied, i.e. rix has eigen- 
value or 1. In the Pauli matrix representation 
this corresponds to eigenstates |0) and |1) of af 
with af\Q) = -|0) and af\l) = |1>. The transfer 
matrix is a product of factors 

eG 

exp{-ehx,i) = exp( — ) 

/exp(-^) \ 

cosh I S sinh | 
^ S sinh I coshf 

V exp(-^)/ 

(14) 

which is a 4x 4 matrix in the Fock space basis |00) , 
|01), 1 10) and |11) of two sites x and x + i. Here 
S — rixAC!fj^iCrfj^2---'^m-i includes the local sign 
r]x,i as well as a non-local string of Pauli matrices 
running over consecutive labels between I and to, 
where I labels the lattice point x and to labels 
X -\- i. Note that S is diagonal in the occupation 
number basis. 

The partition function is now expressed as a 
path integral Z — Sign[n] exp(— S'[n]) over 
configurations of occupation numbers n{x,t) = 
0, 1 on a (3 -I- l)-d space-time lattice of points 
{x,t). The Boltzmann factor is a product of 
space-time plaquette contributions with an action 
S[n{x, t),n{x + i, t), n(x, t-f 1), n{x + i, t+1)] and 
with 

g-s[o,o,o,o] ^ g-s[i,i4,i] ^ e"'^/^ 
g-s[o,i,o,i] =e-s[i,o,i.o] ^^ogj^e 



g-s[o,i,i,o] ^g-s[i,o,o.i] (15) 

All the other Boltzmann weights are zero. 

The sign of a configuration, Sign[n], also is 
a product of space-time plaquette contributions 
Sign[n(x, t), nix + i,t), n{x, t + \), n{x + i,t + 1)] 
with 

Sign[0,0,0,0]) = Sign[0, 1,0,1]) = 
Sign[l,0,l,0]) = Sign[l,l,l,l]) = l, 
Sign[0,l,l,0]) = Sign[l,0,0,l]) = S. (16) 

It should be noted that E gets contributions from 
all lattice points with labels between I and to. 
This seems to make an evaluation of the fermion 
sign rather tedious. Also, it is not a priori obvi- 
ous that Sign[n] is independent of the arbitrarily 
chosen order of the lattice points. Fortunately, 
there is a simple way to compute Sign[n], which 
is directly related to the Pauli exclusion princi- 
ple and is manifestly order-independent. In fact, 
Sign[rt] has a topological meaning. The occu- 
pied lattice sites define fermion world-lines which 
are closed around the Euclidean time direction. 
Of course, during their Euclidean time evolution 
fermions can interchange their positions and the 
fermion world-lines define a permutation of par- 
ticles. The Pauli exclusion principle dictates that 
the fermion sign is just the sign of that permuta- 
tion. If we work with antiperiodic spatial bound- 
ary conditions, Sign[?i] receives an extra minus- 
sign for every fermion world-line that crosses a 
spatial boundary. Figure 1 shows two configura- 
tions of fermion occupation numbers in (1 -I- 1) 
dimensions. The first configuration corresponds 
to two fermions at rest and has Sign[n] = 1. 
In the second configuration two fermions inter- 
change their positions with one fermion step- 
ping over the spatial boundary. If one uses pe- 
riodic spatial boundary conditions this configu- 
ration has Sign[n] = —1. Note that the same 
configuration would have Sign[n] — 1 when an- 
tiperiodic boundary conditions are used. 

Quantities of physical interest are the chiral 
condensate 
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Figure 1. Two configurations of fermion occupa- 
tion numbers m (1 + 1) dimensions. The shaded 
plaquettes carry the interaction. The dots repre- 
sent occupied sites. In the second configuration, 
two fermions interchange their positions. The 
fat line represents a meron-cluster. Flipping this 
cluster changes one configuration into the other 
and changes the fermion sign. 



and the corresponding chiral susceptibility 

x = ^((W>/. 



(18) 



Up to now we have derived a path integral rep- 
resentation for the fermion system in terms of 
bosonic occupation numbers and a fermion sign 
factor that encodes Fermi statistics. The sys- 
tem without the sign factor is bosonic and is 
characterized by the positive Boltzmann factor 
exp(— ^[n]). Here the bosonic model is a quan- 
tum spin system with the Hamiltonian 



H 



q2 q1 



GS S 



3 c3 



x-\-i 



(19) 



where iSJ, = ^cr\ is a spin 1/2 operator associ- 
ated with the lattice site x that was labeled by 
I. From here on, we restrict ourselves to G = 1, 
which corresponds to the antiferromagnetic quan- 
tum Heisenberg model. In the language of the 



spin model, the chiral condensate turns into the 
staggered magnetization 



(20) 



x,t 



The meron-cluster fermion algorithm is based on 
a cluster algorithm for the corresponding bosonic 
model without the sign factor. Bosonic quan- 
tum spin systems can be simulated very efficiently 
with cluster algorithms . The first cluster al- 
gorithm for lattice fermions was described in [^. 
These algorithms can be implemented directly in 
the Euclidean time contimmm i.e. the Suzuki- 
Trotter discretization is not even necessary. The 
decomposition of the lattice into clusters results 
from connecting neighboring sites on each indi- 
vidual space-time interaction plaquette following 
probabilistic cluster rules. A sequence of con- 
nected sites defines a cluster. In this case the clus- 
ters are sets of closed loops. The cluster rules are 
constructed to obey detailed balance. To show 
this constraint we write the plaquette Boltzmann 
factors as 



^-S[n(x,t),n{x+i,t),n(x,t+l),n(x+~i,t+l)] _ 

ASn(x,t),n(x,t+l)^n{x+i,t)M(x+l,t+l) + 
^^n(x,t),l-n(x+l,t)^n{x,t+l)S~-n(x+i,t+l)- 



(21) 



The 5-functions specify which sites are connected 
and thus belong to the same cluster. The quanti- 
ties A and B determine the relative probabilities 
for different cluster break-ups of an interaction 
plaquette. Inserting the expressions from eq.(|TF 
one finds 



A 



'/^,B = sinh- 



(22) 



For plaquette 
configurations [0,0,0,0] or [1,1,1,1] one always 
chooses time-like connections between sites, and 
for configurations [0,1,1,0] or [1,0,0,1] one al- 
ways chooses space-like connections. For configu- 
rations [0,1,0,1] or [1,0,1,0] one chooses time- 
like connections with probability p — A/{A 
-B) = 2/[l -I- exp(e/2)] and space-like connections 
with probability 1 — p — B/{A -\- B). Indeed, 
this is the algorithm that was used in [|o|. It 
is extremely efficient, has almost no detectable 
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autocorrelations, and has a dynamical exponent 
for critical slowing down that is compatible with 
zero. The cluster rules are illustrated in table 1. 
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Table 1 



Cluster break-ups of various plaquette configura- 
tions together with their relative probabilities A 
and B. The dots represent occupied sites and the 
fat lines are the cluster connections. 



Eq. (gl|) can be viewed as a representation of the 
original model as a random cluster model. The 
cluster algorithm operates in two steps. First, a 
cluster break-up is chosen for each space-time in- 
teraction plaquette according to the above prob- 
abilities. This step effectively replaces the orig- 
inal Boltzmann weight of a plaquette configura- 
tion with a set of constraints represented by the 
(5-functions associated with the chosen break-up. 
The constraints imply that occupation numbers 
of connected sites can only be changed together. 
Second, every cluster is flipped with probability 



1/2. When a cluster is flipped the occupation 
numbers of all sites that belong to the cluster are 
changed. Eq.(p^) ensures that the cluster algo- 
rithm obeys detailed balance. 

Let us now consider the effect of a cluster flip 
on the fermion sign. Each cluster can be charac- 
terized by its effect on the fermion sign indepen- 
dent of the orientation of all the other clusters. 
We refer to clusters whose flip changes Sign[n] as 
merons, while clusters whose flip leaves Sign[n] 
unchanged are called non-merons. The flip of a 
meron-cluster permutes the fermions and changes 
the topology of the fermion world-lines. The 
number of merons in a configuration is always 
even. An example of a meron-cluster is given in 
figure 1. 

The meron concept alone allows us to gain an 
exponential factor in statistics. Since all clusters 
can be flipped independently, one can construct 
an improved estimator for Sign[n] by averaging 
analytically over the 2^<^ configurations obtained 
by flipping the Nq clusters in the configuration 
in all possible ways. For configurations that con- 
tain merons the average Sign[n] is zero because 
flipping a single meron leads to a cancellation of 
contributions of ±1. Hence, only the configura- 
tions without merons contribute to (Sign). The 
vast majority of configurations contains merons 
and now contributes an exact to (Sign) instead 
of a statistical average of contributions ±1. In 
fact, one can show that the contributions from the 
zero-meron sector are always positive, such that 
Sign[n] is for configurations containing meron- 
clusters and 1 in the zero-meron sector. Accord- 
ing to the previous discussion, this method solves 
one half of the fermion sign problem. Before we 
can solve the other half of the problem we must 
discuss improved estimators for the physical ob- 
servables. 

Let us consider an improved estimator for 
(^'5'[n])^Sign[n] which is needed to determine the 
chiral susceptibility x- The total chiral conden- 
sate, ^'^'[ti] — J2c is a sum of cluster con- 
tributions 

= ^ E (23) 
When a cluster is flipped, its condensate contri- 
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bution changes sign. In a configuration without 
merons, where Sign[n] = 1 for all relative cluster 
flips, the average of (^"I'[n])^Sign[n] over all 2^"^ 
configurations is J^c l^^cP- For configurations 
with two merons the average is 2|5'5'c'i I |^^C2 1 
where Ci and C2 are the two meron-clusters. 
Configurations with more than two merons do not 
contribute to (^'5'[n])^Sign[?i]. The improved es- 
timator for the susceptibility is hence given by 

^ (Ec l^^cP^jv,o + 2\mc, I I^^C, \Sn,2) 

(24) 

where N is the number of meron-clusters in a con- 
figuration. Thus, to determine x oii^ must only 
sample the zero- and two-meron sectors. 

The probability to find a configuration with- 
out merons is exponentially small in the space- 
time volume since it is equal to (Sign). Thus, 
although we have increased the statistics tremen- 
dously with the improved estimators, without a 
second step one would still need exponentially 
large statistics to accurately determine x- For- 
tunately, the numerator in equation (p^ ) receives 
contributions from the zero- and two-meron sec- 
tors only, while the denominator gets contribu- 
tions only from the zero-meron sector. One can 
hence restrict oneself to the zero- and two-meron 
sectors and never generate configurations with 
more than two merons. This enhances both the 
numerator and the denominator by a factor that 
is exponentially large in the volume, but leaves 
the ratio of the two invariant. One purpose of 
the second step of the meron-cluster algorithm is 
to eliminate all configurations with more than two 
merons. To achieve this, we start with an initial 
configuration with zero or two merons. For ex- 
ample, a completely occupied configuration has 
no merons. We then visit all plaquette interac- 
tions one after the other and choose new pair 
connections between the four sites according to 
the above cluster rules. If the new connection in- 
creases the number of merons beyond two, it is 
not accepted and the old connection is kept for 
that plaquette. This procedure obeys detailed 
balance because configurations with more than 
two merons do not contribute to the observable 
we consider. This simple reject step eliminates 



almost all configurations with weight and is the 
essential step to solving the second half of the 
fermion sign problem. 

We have simulated the staggered fermion 
model with G — 1 on antiperiodic spatial vol- 
umes with L — 4:,6, 16 and at various in- 
verse temperatures /? G [0.5, 1.2] which includes 
the critical point. In the Euclidean time direction 
we have used M = 4, i.e. 24 time-slices. In all 
cases, we have performed at least 1000 thermal- 
ization sweeps followed by 10000 measurements. 
One sweep consists of a new choice of the cluster 
connections on each interaction plaquette and a 
flip of each cluster with probability 1/2. 
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Figure 2. The probability of having a certain num- 
ber of merons for spatial size L ~ A and 10 at 
(3 = 0.948. 
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Figure 3 shows the probabihty to have a cer- 
tain number of merons in an algorithm that sam- 
ples all meron-sectors. For small volumes the 
zero-meron sector and hence (Sign) is relatively 
large, while multi-meron configurations are rare. 
On the other hand, in larger volumes the vast 
majority of configurations has a large number of 
merons and hence (Sign) is exponentially small. 
For example, an extrapolation from smaller vol- 
umes gives a rough estimate (Sign) « 10~^° on 
the 16'^ lattice at /? = 0.948. Hence, to achieve 
a similar accuracy without the meron-cluster al- 
gorithm one would have to increase the statistics 
by a factor 10''°, which is obviously impossible in 
practice. 

To study the critical behavior in detail, we have 
performed a finite-size scaling analysis for x fo- 
cusing on a narrow range (3 € [0.9, 0.98] around 
the critical point. Since a Z(2) chiral symmetry 
gets spontaneously broken at finite temperature 
in this (3 -I- l)-d model, one expects to find the 
critical behavior of the 3-d Ising model. For the 
3-d Ising model the critical exponents are given 
by 1/ = 0.630(1) and -(/v = 1.963(3). Fitting to 
our data, we find v = 0.63(4) and j/v = 1.98(2), 
which indicates that the chiral transition of the 
staggered fermion model is indeed in the 3-d Ising 
universality class. In figure 3 we have taken the 
values of the critical exponents from the 3-d Ising 
model and we have plotted xl^'^^^ ^ function 
oiy = {P — l3c)L^^'^ ■ Indeed all data collapse onto 
one universal curve. 

4. Quantum Antiferromagnets in a Mag- 
netic Field"' 

As another problem of physical interest we con- 
sider antiferromagnetic quantum spin ladders in 
an external uniform magnetic field that competes 
with the spin-spin interaction. The competing 
interactions lead to a severe sign problem in nu- 
merical simulations. The most efficient algorithm 
for simulating quantum spin systems in the ab- 
sence of an external magnetic field is the loop 
cluster algorithm |^,|l^,^ which can be imple- 
mented directly in the Euclidean time continuum 
0. This algorithm also underlies our staggered 
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Figure 3. Finite-size scaling behavior of the chiral 
susceptibility x- The data for various spatial sizes 
L = 4, 6, 8 and 10 fall on one universal curve. 



fermion simulations and was described in detail 
in section 3. In the presence of a magnetic field 
pointing along the quantization axis of the spins, 
the Z(2) symmetry that allows the clusters to be 
flipped with probability 1/2 is explicitly broken. 
As a consequence, for large values of the magnetic 
fleld some clusters can only be flipped with very 
small probability and the algorithm becomes in- 
efficient. An alternative strategy is to choose the 
spin quantization axis perpendicular to the direc- 
tion of the magnetic field. Then all clusters can 
still be flipped with probability 1/2. This indeed 
leads to an efficient algorithm for ferromagnets 
in an external uniform magnetic field. Unfortu- 
nately, for antiferromagnets i.e. when the mag- 
netic field competes with the spin-spin interaction 
— this formulation of the problem leads to a very 
severe sign problem. Here we show how this sign 
problem can be solved completely using a meron- 
cluster algorithm. 

Antiferromagnetic spin ladders — sets of sev- 
eral transversely coupled quantum spin chains — 
are interesting condensed matter systems which 
interpolate between single 1-d spin chains and 
2-d quantum antiferromagnets. The ladders are 
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spatially quasi 1-d systems whose low-energy dy- 
namics are governed by (1 -I- l)-d quantum field 
theories. We consider spin ladders in an external 
uniform magnetic field i?, which corresponds to 
a chemical potential fj, = B/c in the correspond- 
ing (1 -|- l)-d quantum field theory. Chakravarty, 
Halperin and Nelson used a (2 -I- l)-d effective 
field theory to describe the low-energy dynamics 
of spatially 2-d quantum antiferromagnets JlH . 
Chakravarty has applied this theory to quantum 
spin ladders with a sufficiently large even num- 
ber of coupled spin 1/2 chains [|l2|. We consider 
spin ladders with the same value of the antiferro- 
magnetic coupling along and between the chains. 
These systems are described by the action 

5[e| = dt dx dy ^[d^S-dxe 
Jo Jo Jo ^ 

+ dyg- dye + -^dte- dte\. (25) 

Here e{x,y,t) is a unit-vector field, ps is the 
spin stiffness and c is the spin-wave velocity. 
The coupled spin chains are oriented in the spa- 
tial a;-direction with a large extent L, while the 
transverse y-direction has a much smaller extent 
L' <ti L. Here we consider spin ladders with peri- 
odic boundary conditions in the transverse direc- 
tion, and we limit ourselves to an even number of 
coupled spin 1/2 chains. The effective action for 
a ladder with an odd number of coupled chains 
would contain an additional topological term. 

When the spin ladder is placed in a uniform 
external magnetic field B, the field couples to a 
conserved quantity — the uniform magnetization. 
Hence, on the level of the effective theory, the 
magnetic field plays the role of a chemical poten- 
tial — i.e. it appears as the time-component of 
an imaginary constant vector potential. The or- 
dinary derivative 9te is replaced by the covariant 
derivative dte + iB x e. For a sufficiently large 
even number of coupled chains (L' ^ c/ ps) the 
ladder system undergoes dimensional reduction 
to the (1 + l)-d 0(3) symmetric quantum field 
theory with the action 

S[e\^ dt dx t!±-[d^e-dxe 
Jo Jo ^ 



+ ^(dte + iB X e) ■ (dte + iB X e)]. (26) 

The effective coupling constant is given by 1 /g'^ = 
PsL' /c and the magnetic field appears as a chem- 
ical potential of magnitude p, — B/c. 

We consider a system of quantum spins 1/2 
on a d-dimensional cubic lattice with site label 
X and with periodic spatial boundary conditions. 
In particular, we are interested in ladder systems 
on a 2-d rectangular lattice of size L x L' with 
L L' . The spins located at the sites x are 
described by operators with the usual com- 
mutation relations 

[Sl,Sl]=i6xyO.,kSt (27) 
The Hamilton operator 

H = jJ2Sx-S,^^i-B-J2Sx, (28) 

x.i X 

with J > 0, couples the spins at the lattice sites 
X and X -f i, where j is a unit- vector in the i- 
direction. A path integral representation of the 
partition function can be derived in complete 
analogy to the staggered fermion model of section 
3. In this case, one sums over configurations of 
spins s(x,t) =t, i on a (d-l- l)-dimensional space- 
time lattice of points {x,t). Again, the Boltz- 
mann factor is a product of space-time plaquette 
contributions with Sign[s(a;, t), s{x + i, t), s{x, t + 
1), s{x + i,t+ 1)] given by 

Sign[T,T,T,T] = Sign[t,i,T,i] = 
Sign[i,t,i,T] = Sign[i,i,i,i] = l 

Sign[t,i,i,T] = Sign[i,t,Ta] = -l. (29) 

When the magnetic field points perpendicular to 
the spin quantization axis, time-like bond contri- 
butions to the Boltzmann factor also arise which 
have 

g-s[T,Tl ^ e-su,il ^ cosh(eB/2) 

e-^[T,i] ^ g-s[i,T] ^ sinh(e5/2). (30) 

Figure 4 shows two spin configurations in (1 + 1) 
dimensions. The first configuration is completely 
antiferromagnetically ordered and has Sign[s] = 
1. The second configuration contains one in- 
teraction plaquette with configuration [i,!:?;-!] 
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which contributes Sign[|, t, t, J,] = — 1. In ad- 
dition, there are two time-hke interaction bonds 
with configurations [|,t] and [t,i]- 
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Figure 4. Two spin configurations in (1 + 1) 
dimensions. The shaded plaquettes and shaded 
time-like bonds carry the interaction. Filled dots 
represent spin up and open circles represent spin 
down. The second configuration has Sign[s\ — 
— 1. A meron- cluster is represented by the fat 
black line. Flipping this cluster changes one con- 
figuration into the other and changes Sign[s] . 



gered fermions. Just hke the space-time plaque- 
tte terms, the time-hke bond Boltzmann factors 
are expressed as 



-S[s(x,t),s(x,t+l)] 



C5. 



s{x .t) ^s{x .t^l) 



D. 



(32) 



The probabihty to connect spins with their time- 
hke neighbors is C /{C + D). The spins remain 
disconnected with probabihty F) /(C -f D). In- 
serting the expressions from eg. (pOj) one obtains 



C + D = cosh(eS/2), C = sinh(eB/2). 



(33) 



The cluster rules for the space-time plaquette 
terms are the same as in table 1. The rules for 
the time-like bonds are illustrated in table 2. 



configuration 


break-ups 
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D 



Table 2 

Cluster break-ups of time-like bond configurations 
together with their probabilities C and D. Filled 
dots represent spin up, open circles represent spin 
down, and the fat black line is the cluster connec- 
tion. 



The central observable of our study is the uni- 
form magnetization 

M = Y.^- (31) 

The expectation value of the magnetization {M) 
in the direction of the magnetic field is non-zero. 

The cluster algorithm for the quantum anti- 
ferromagnet is very similar to the one for stag- 



The above cluster rules were first used in a 
simulation of the Heisenberg antiferromagnet ||l^ 
in the absence of a magnetic field. In that case 
there is no sign problem. Then the correspond- 
ing loop-cluster algorithm is extremely efficient 
and has almost no detectable autocorrelations. 
When a magnetic field is switched on the situ- 
ation changes. When the magnetic field points 
in the direction of the spin quantization axis (the 
3-direction in our case) there is no sign problem. 
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However, the magnetic field then explicitly breaks 
the Z(2) flip symmetry on which the cluster al- 
gorithm is based, and the clusters can no longer 
be flipped with probability 1/2. Instead the flip 
probability is determined by the value of the mag- 
netic field and by the magnetization of the clus- 
ter. When the field is strong, flips of magnetized 
clusters are rarely possible and the algorithm be- 
comes very inefficient. To avoid this, we have 
chosen the magnetic field to point perpendicular 
to the spin quantization axis. In that case, the 
cluster flip symmetry is not affected by the mag- 
netic field, and the clusters can still be flipped 
with probability 1/2. However, now a severe sign 
problem arises. 

One can construct an improved estimator for 
the magnetization (M) , which gets non- vanishing 
contributions only from the zero-meron sector. 
Hence, it is unnecessary to generate any configu- 
rations that contain meron-clusters. In practice 
it is advantageous to occasionally generate con- 
figurations containing merons even though they 
do not contribute to our observable, because this 
reduces the autocorrelation times. 



We have performed numerical simulations with 
the meron-cluster algorithm for various quantum 
antiferromagnets in a uniform magnetic field B. 
To demonstrate the efficiency of the algorithm, we 
have compared it with the standard loop-cluster 
algorithm. In case of the loop algorithm the mag- 
netic field points in the direction of the spin quan- 
tization axis. In the meron-cluster algorithm, on 
the other hand, the magnetic field is perpendic- 
ular to the spin quantization axis and all clus- 
ters can still be flipped with probability 1/2. Of 
course the sign problem arises, but it is solved 
completely by the meron-cluster algorithm. Fig- 
ure 5 compares the thermalization behavior of the 
magnetization of a 2-d Heisenberg antiferromag- 
net on an 8 X 8 lattice at /3J = 10 with M = 100 
for the loop-cluster algorithm and the meron- 
cluster algorithm. At B = J the loop-cluster 
algorithm needs more than 100000 equilibration 
sweeps, while the meron-cluster algorithm has no 
thermalization problem. 
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Figure 5. Thermalization of the magnetization 
for the loop- cluster algorithm versus the meron- 
cluster algorithm. For B = J the loop algorithm 
takes more than 100000 sweeps to reach equilib- 
rium while the meron-cluster algorithm has no 
thermalization problem. 
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Figure 6. Magnetization density as a function of 
the magnetic field. 



Figure 6 shows the magnetization density of an- 
tiferromagnetic quantum spin ladders consisting 
of four coupled spin 1/2 chains (i.e. L' /a = 4). 
The cases L = 20, /3J = 15, M = 200 and 
L = 40, PJ = 24, M = 300 have been considered. 
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In the infinite volume, zero temperature limit one 
expects the magnetization to vanish for B < B^- 
The critical magnetic field Be corresponds to the 
critical chemical potential n = B^jc = m, deter- 
mined by the mass gap m of the (1 + l)-d 0(3) 
model. The expected behavior is indeed consis- 
tent with the numerical data. 

5. Strongly correlated Electron Systems^ 

Quantum Antifcrromagncts may turn into 
high-temperature superconductors when they are 
doped with additional charge carriers. These sys- 
tems are believed to be well described by the Hub- 
bard model Hamiltonian 

x,i,s='\,l 

X X 

(34) 

Here s "^x.s arc creation and annihilation 
operators for electrons with spin s =t, i hopping 
on the sites a; of a quadratic lattice. The chemical 

potential /j, coiiplcs to the occupation number per 

site rix = nx,'\ + n^a = *x,t*x,T + *x,i*^.i ^"^"^ 
is used to dope the system away from half-filling. 
So far, high-temperature sup{^rc;onductivity has 
not been demonstrated convincingly in numeri- 
cal simulations due to a very severe fermion sign 
problem. It is natural to ask if a mcron-cluster 
algorithm can be used to solve this problem. Un- 
fortunately, as it stands the meron concept does 
not apply directly to the Hubbard model because 
the clusters influence each other in their effect on 
the fermion sign. 

Meron-cluster algorithms can solve the sign 
problem when two conditions are satisfied. First, 
the clusters must be independent in their effect on 
the sign in order to allow the construction of im- 
proved estimators. Second, it must always be pos- 
sible to flip the clusters into a reference configura- 
tion with a positive sign. In case of the Hubbard 
model a completely antiferromagnetically ordered 
reference configuration suggests itself. In order 
to construct an efficient meron-cluster algorithm 

^Based on a talk presented by J. Cox 



based on that reference configuration, one must 
modify the original Hubbard model Hamiltonian. 
Of course, one would like to maintain the symme- 
try properties of the original Hamiltonian in order 
to stay in the same universality class. An obvious 
symmetry is the SU (2)^ spin rotational symmetry 

that is generated hy S = \ Y^x,s'^t,s^s,s''^ x,s' ■ 
Another less obvious 5f7(2)„ symmetry is gener- 
ated by N with 

7V+ = ^(-lf i+^^*+^*+^,Ar- = {N+)+, 

X 

^' = ^E«T*-.T + *+i*.,i)- (35) 

X 

We have systematically investigated the space 
of SU{2)s (8) SU{2)n nearest neighbor interaction 
Hamiltonians for which an efficient meron-cluster 
algorithm can be constructed. One example of 
such a Hamiltonian is 

x,i,s=1,l 

X (n, + n^+j - 2)2 

x,i 

-|-C/^(n^,t - ^){n^,i - i) 

X 

+ y^x,i[(n.-l)2 + K+j-l)2 

- 2(n, - l)2(n,+. - If] - mE("- - !)• (36) 

X 

In particular, an additional antiferromagnetic 
coupling JSx ■ S^_y_i arises, which is constrained 
by J > 2t. Furthermore, t' = t, U > J 
and U > 8V. Indeed the above system can be 
simulated reliably with a very efficient meron- 
cluster algorithm. Figure 7 shows the density of 
charge carriers away from half-filling as a func- 
tion of the chemical potential for various temper- 
atures. At /i « 0.45 there is a strong first order 
phase transition at which the system undergoes 
phase separation. So far we have not found high- 
temperature superconductivity. This result may 
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be a consequence of the additional antiferromag- 
netic coupling. However, it could also indicate 
that even the original Hubbard model does not 
display high-temperature superconductivity. 




O:/;j=1.0 
□:/SJ = E.O 
*:/JJ = 3.0 
*:(SJ = 4.0 



Figure 7. Particle density as a function of chem- 
ical potential for four different temperatures. 



6. The 2-d 0(3) Model at non-zero Chem- 
ical Potential* 

Another severe sign problem arises in simula- 
tions of field theories at non-zero chemical po- 
tential. Obviously, understanding QCD at non- 
zero baryon density is of great experimental and 
theoretical importance. Here we study a simpler 
toy model for QCD — the 2-d 0(3) model. Like 
QCD, this model is asymptotically free and has a 
non-perturbatively generated mass gap, instan- 
tons and 9 vacua. In addition, it also suffers 
from a very severe sign problem when coupled 
to a chemical potential. As already discussed in 
section 4, the action of the 2-d 0(3) model with 
chemical potential jl takes the form 




+ {d2e + iflxe)-{d2e + iflxe)]. (37) 

Note that in this case the chemical potential is a 
vector because it couples to a non-Abelian con- 
served charge. 

*Based on two talks given by B. Scarlet and U.-J. Wiese 



Due to the use of discrete variables, a meron- 
cluster algorithm can be applied to the D-theory 
formulation of the 2-d 0(3) model. Instead of 
performing a path integral over classical fields, 
in the D-theory formulation of field theory fl^ 
collective excitations of discrete quantum vari- 
ables which undergo dimensional reduction play 
the role of an effective classical field. For exam- 
ple, in D-theory the 2-d 0(3) model arises from 
dimensional reduction of the (2 -I- l)-d antifer- 
romagnetic quantum Heisenberg model with the 
Hamiltonian 

H = jY,S.-S..+1- (38) 

x.i 

At zero temperature this model develops a stag- 
gered magnetization and the SO{3) spin rota- 
tional symmetry breaks spontaneously down to 
5*0(2). The low-energy dynamics of the corre- 
sponding Goldstone bosons — the magnons or 
spin waves — are described by the chiral pertur- 
bation theory action 

S[e\^ j^\t j d^x^[die-die+^dte-dt^.[m) 

Here ps is the spin stiffness and c is the spin wave 
velocity. At low temperatures the correlation 
length of the system diverges as ^ cx exp(27rps/3). 
Consequently, the extent (3 of the Euclidean time 
direction becomes negligible compared to ^ and 
the system undergoes dimensional reduction. The 
dimensionally reduced theory is the 2-d 0(3) 
model 

S[e\ = j d^x^d,e-d,e (40) 

without chemical potential. In order to incorpo- 
rate a non-zero chemical potential, the D-theory 
Hamiltonian must be modified to 

X 

+ ^(^^ ^;+2^' + s^K^2^'') + ^'^'+2]- (41) 

This system dimensionally reduces to the 2-d 
0(3) model with chemical potential of eq.(|37|). 

The above Hamiltonian can be simulated with 
a meron-cluster algorithm very much like the one 
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for quantum antiferromagnets in a magnetic field. 
The resulting particle density as a function of 
the chemical potential is shown in figure 8. The 
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Figure 8. The particle density as a Junction of the 
chemical potential in the 2-d 0(3) model. The 
Monte Carlo data are consistent with analytic 
predictions at low density and at infinite volume 
represented by the dashed lines. 



Monte Carlo results show the expected onset be- 
havior when the chemical potential exceeds the 
mass gap m — 1/^. In the zero temperature limit 
no particles would be produced for fx < m. At 
finite temperatures, on the other hand, a small 
particle density which can be computed analyt- 
ically exists in this region. As shown in figure 
8, the Monte Carlo data are consistent with this 
analytic prediction. For larger values of finite 
temperature effects are small and one can com- 
pare with the thermodynamic Bethe ansatz solu- 
tion of the 2-d 0(3) model, which again is consis- 
tent with the numerical data. 

The steps taken here for the 2-d 0(3) model can 
potentially all be generalized to QCD. In the D- 
theory formulation, QCD is described as a quan- 
tum link model , and a chemical potential can 
be included exactly as in eq.(|4l|). The quarks 
in quantum link QCD appear as domain wall 
fermions. Hence, to simulate QCD with chemical 
potential along these lines, one must first con- 
struct a meron-cluster algorithm for domain wall 
fermions. 



7. Domain Wall Fermions ^ 

For simplicity, let us consider the Hamiltonian 
of free (2 -j- l)-d domain wall fermions, which de- 
scribes 2-d chiral fermions bound to the walls 

x,i 

+ M^vE't^3*x. (42) 

X 

Here and are quark creation and annihi- 
lation operators obeying canonical anticommuta- 
tion relations. Following ||l^, the partition func- 
tion is written as 

Z= (0|exp(-/3i7)|0), (43) 

where |0) is a particular fermion Fock state. Tak- 
ing the expectation value in that state implies 
that there are no left-handed quarks at a;3 = 0, 
and no right-handed quarks at X3 = (3. 

Like for the Hubbard model, we have thus far 
not found an efficient meron-cluster algorithm for 
the original domain wall fermion Hamiltonian. 
Instead, we ask for which Hamiltonians an effi- 
cient algorithm can be constructed. In order to 
stay in the same universality class, we demand 
that these Hamiltonians have the same symme- 
try properties as the original domain wall fermion 
Hamiltonian. The relevant symmetries are charge 
conjugation C, parity P, and lattice rotations R 
by 90 degrees. As in the staggered fermion case 
of section 3, we decompose the Hamiltonian as 
H = j hx,i and we consider the transfer ma- 
trix Ti — exp(— eft.a;_,;). In the Hilbert space, the 
symmetries C, P and R are represented by unitary 
transformations Uc^ Up and Ur. The symmetry 
requirements on the (2 + l)-dimensional transfer 
matrix thus take the form 

= U+pT,Uc^p, T2 = U+TiUr. (44) 

In order to stay in the right universality class, it 
is equally important not to have any additional 

^Based on a talk presented by C. Gattringer 
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symmetries that are broken in the original do- 
main wall fermion Hamiltonian. This is partic- 
ularly important for removing doubler fermions. 
Like for the Hubbard model, we have system- 
atically investigated all symmetry requirements, 
and we have identified the most general near- 
est neighbor interaction Hamiltonian for which 
an efficient meron-cluster algorithm can be con- 
structed. Compared to previous cases an addi- 
tional complication arises. The 7-matrices give 
rise to contributions ±i to the total Boltzmann 
weight. Hence, the flip of a general cluster may 
not only change the sign — it may change the 
complex phase of the Boltzmann weight by a fac- 
tor i. Consequently, besides the sign-changing 
meron-clusters there are i-on clusters whose flip 
changes the phase of the Boltzmann weight by a 
factor i. 

At present, we have not fully explored the 
space of efficient cluster algorithms for domain 
wall fermions. The necessary modifications of the 
standard Hamiltonian are similar to the ones re- 
quired in the Hubbard model. It remains to be 
seen if chiral fermions bound to the walls appear 
in the modified models that can be simulated with 
efficient meron-cluster algorithms. 

8. Conclusion 

In conclusion, the meron concept provides us 
with a powerful algorithmic tool — the meron- 
cluster algorithm — which can lead to a com- 
plete solution of severe sign problems. Here we 
have demonstrated that this algorithm allows us 
to simulate staggered fermions with an unusu- 
ally small number of flavors, certain strongly cor- 
related electron systems, quantum antiferromag- 
nets in a magnetic field, as well as the 0(3) model 
at non-zero chemical potential. The next chal- 
lenge is to construct meron-cluster algorithms for 
QCD at non-zero baryon density and for systems 
that show high-temperature superconductivity. 
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